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We use a quantum Monte Carlo method to study the ground state and thermodynamic phase transitions of 
the spin supersolid phase in the S = 1 Heisenberg model with uniaxial anisotropy. The thermal melting of 
the supersolid phase shows unqiue signatures in experimentally measurable observables. This Hamiltonian is 
a particular case of a more general and ubiquitous model that describes the low energy spectrum of a class of 
isotropic and frustrated spin systems. We also discuss some alternative realizations of spin supersolid states in 
real magnets. 
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The supersolid (SS) state of matter has attracted great in- 
terest lately following the experiments of Kim and Chan on 
solid 4 He. While it is still unclear whether it can be stabilized 
in the continuum, there are several numerical studies which 
show that a SS phase can be realized in the presence of a peri- 
odic potential or underlying lattice for bosons QSHl as well 
as spins lit]. The SS state is easier to stabilize on a lattice 
because the lattice parameter of the "solid phase" or charge 
density wave cannot relax to any arbitrary value (it has to be 
an integer multiple of the underlying lattice parameter). In this 
work we have studied a class of spin-SS on cubic lattices, fo- 
cusing primarily on the unique signatures of the thermal melt- 
ing of the SS in experimentally measurable observables. We 
also discuss different conditions under which a spin-SS can be 
realized in real spin compounds. 

The minimal spin model that has a thermodynamically sta- 
ble supersolid ground state on a bipartite lattice is the S = 
1 Heisenberg model with uniaxial single-ion and exchange 
anisotropics and an external magnetic field: 

h h = jJ2i s i s ? + s i s ! + AS i S D +J2( DS i 2 - BS ^ 

<ij> i 

(1) 

where (i, j) indicates that i and j are nearest neighbor sites, 
D is the amplitude of the single ion-anisotropy and A de- 
termines the magnitude of the exchange uniaxial anisotropy. 
Note that although the exchange interaction is anisotropic, the 
longitudinal ( J) and transverse (A) couplings are both AFM 
(positive). Henceforth, J is set to unity and all the parameters 
are expressed in units of J. 

The ground state properties of the above model on a square 
lattice were studied in detail previously |@]. For D, A > 1, the 
ground state is supersolid over a finite range of applied field 
B. In this work, we report the ground state and thermody- 
namic properties of (HJ on a cubic lattice. While the quantum 
phase diagram remains qualitatively unchanged, the thermo- 
dynamic properties are different in three dimensions (3Ds) be- 
cause the condensate (XY ordered antiferromagnetic phase) 
extends to finite temperatures. Additionally, the melting of 
this phase belongs to the XY universality class as opposed to 
Kosterlitz-Thouless (KT) type in two dimensions (2Ds). This 
has important consequences in any putative experimental de- 



tection of the SS phase, as we shall discuss below. 

The Stochastic Series expansion (SSE) quantum Monte 
Carlo (QMC) methodQ is used to study the Hamiltonian (Q} 
on cubic lattices N = L x L x L/2, with 8 < L < 16. To 
characterize the different phases, we compute the longitudinal 
component of the staggered static structure factor (SSSF), 

5Z2 (Q) = ^E e " q ' (r3 ~ rfc) ^>' Q = (tTj 7r), 

(2) 

and the spin stiffness, p s , defined as the response of the system 
to a twist in the boundary conditions. 5 Z2 (Q) measures the 
extent of diagonal (Ising like) long-range order (LRO) at the 
ordering wave vector Q = (it, ir, ir), while the stiffness (su- 
perfluid density in particle language 18D), indicates the pres- 
ence of XY (off-diagonal) LRO (this is not true for D< 3). 
In 3D, the superfulid density is identical to the condensate 
fraction. In the simulations, the stiffness is obtained from 
the winding numbers of the world lines along the three axes: 
p s = (W* + W* + W 2 2 )/3/3i. 

Ground state (GS) phases As the field, B, is varied, the GS 
of (Q3 goes through a succession of phases, including spin- 
gapped Ising-ordered (IS) and gapless XF-ordered (XY) 
states. The IS phase is marked by a diverging value of 
S ZZ (Q) oc N in the thermodynamic limit, whereas a finite 
p s characterizes the gapless XY ordered phase. A spin SS 
phase is identified by a finite value of both S ZZ (Q)/N and 
Ps ifioll . Both quantities are always finite for finite size sys- 
tems and estimates for N — > oo are obtained from finite-size 
scaling. 

Fig. [TJ shows the quantum phase diagram as a function of 
magnetic field, B, for D — 3.0 and A = 6.0 - it is qualita- 
tively similar to that obtained in 2Dsj6t]. The m z (B) curve 
features two prominent plateaus corresponding to different IS 
phases. For small B, the GS is a gapped AFM solid (IS 1 ) 
with no net magnetization. The stiffness, p s , vanishes in the 
thermodynamic limit, while S ZZ (Q)/N 1 with the spins 
primarily in the S? = ±1 states in the two sublattices. At a 
critical field, B c \, there is a second order transition to a state 
with a finite fraction of spins in the Sf = state. These 
S z = "particles" Bose-Einstein condense (BEC) to give the 
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FIG. 1 : (Color online) Quantum phase diagram of Hh (Eq. [TJ for 
D — 3.0 and A = 6.0. The upper panel shows the magnetization 
as a function of field B. The SS phase appears between the two 
Ising-like phases IS1 and IS2. At higher fields, there is a first order 
transition to a pure XY-AFM phase. The lower panel shows the 
stiffness and the longitudinal component of the SSF. The SS phase 
has finite values of both observables. 



FIG. 2: (Color online) Two-step melting of the SS for parameters in 
Fig.l and B / J=\6.5 (solid line in Fig.l). The top panel shows the 
vanishing of the two orders at different temperatures. The disappear- 
ance of superfluidity is accompanied by an unusual increase in the 
solid order parameter. The lower panel shows the signatures in the 
specific heat at the two transitions. 



GS a finite stiffness. The diagonal order is reduced but re- 
mains finite as well. The resulting GS thus has simultaneous 
long-range diagonal and off-diagonal order; in other words, it 
is a spin-supersolid. The complete phase diagram consists of 
a second gapped Ising phase (IS2) with diagonal order (all the 
spins in the S z = — 1 sublattice are flipped to S z = 0) and an 
XY phase at very high fields with pure off-diagonal ordering. 

Finite temperature transitions Finite temperature proper- 
ties of the SS has previously been examined for hard core 
bosons [|3l [Tl, 12] and S = \ spins on a bilayer|13|. The 
melting of the SS phase proceeds via two steps - the super- 
fluid order disappears at a lower temperatute whereas the solid 
order persists up to a higher temperature. In 2Ds, the contin- 
uos U(l) symmetry cannot be broken at T > and the SS has 
only a quasi long-range off-diagonal order for T < Tkt- The 
vanishing of the spin stiffness occurs via a KT transition. In 
contrast, true long-range off-diagonal order in the SS persists 
to finite temperatures in 3Ds and the melting of the superfluid 
order belongs to the XY universality class. The solid order 
disappears at a higher temperature via an Ising-like transition. 

The results of simulations of thermal transitions associated 
with the SS phase are shown in fig|2] The top panel shows 
the variation of the solid and superfluid order parameters as 
a function of temperature. At low temperatures, both order 
parameters are finite. With increasing T, the SS "melts" into 
a pure solid. The disappearence of SF order is marked by an 



enhancement in the solid order. This apparently anomalous 
behavior reflects the fact that in the SS phase, the solid order 
is partially suppressed by the co-existing SF order. The lon- 
gitudinal component of the SSF is accessible experimentally 
by neutron scattering and its non-monotonic behavior at the 
onset of superfluid order can serve as an important signature 
of the SS phase. The three dimensionality of the model im- 
plies that the two transitions should be accompanied by spe- 
cific heat anomalies at the corresponding temperatures. The 
XY transition will manifest itself as a A-anomaly while the 
Ising-like solid melting will be marked by a cusp. Indeed we 
find clear signatures of the two transitions in the calculated 
specific heat (lower panel of fig|2]i. While both the peaks are 
rounded by finite-size effects, their positions coincide unam- 
biguously with the melting of the superfluid and Ising orders. 
Since it is one of the most readily measurable observables, 
having clear signatures in the specific heat is of great experi- 
mental relevance. 

Next we discuss the relevance of these results for finding 
a SS phase in real magnets. The magnetic properties of spin 
compunds with spin-orbit interaction much smaller than the 
crystal field splitting can be adequately described by a U(l ) 
invariant model (although this invariance is never perfect) II 1 011 . 
The transition metal magnetic ions belong to this class. On the 
other hand, the exchange anisotropy is typically very small in 
these compunds. The above model with large A is not di- 
rectly applicable to this class of real quantum magnets. We 
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shall show below that under appropriate conditions, an effec- 
tive uniaxial exchange anisotropy can be generated in the low- 
energy subspace of a model with (realistic) isotropic interac- 
tions. To this end we consider coupled layers of dimers with 
only isotropic (Heisenberg) AFM interactions - an intra-dimer 
exchange Jo and weaker inter-dimer /hwfrafe J couplings J\ 
and J 2 (see Figs.[3{a) and[3jb)): 



Hd — Jo Si + • S,_ + Ji S- la ■ Sj Q 

+ J 2 ^2 S ia ■ Sj S - B Sia- 

<ij>>" ia 



(?) 



The index a = ± denotes the two spins on each dimer. 

For 5 = 1 dimers, the low energy subspace of Ho (for 
Ji , J 2 -C Jo) consists of the singlet, the S z = 1 triplet and 
the S z = 2 quintuplet (see Fig. [3ja)). The low energy effec- 
tive model that results from restricting Ho to this subspace 
supports a field-induced supersolid phase on a bipartite lattice 
for J > z(Ji + Ja)/2 and J > z{J\ - J 2 )/2 (z is the 
co-ordination number of the lattice) as was shown in Ref.|01. 



(a) 




1 J 2 

Square dimer lattice 



(c) 

S=l 



operator b\bi and the parameters of the effective model are 
expressed in terms of those of the original Hamiltonian Hrj 
as t={Jx- J 2 )/2, V = (Ji + J 2 )/2 and ti = -J + B. 
On many frustrated lattices, this model contains a SS phase 
in its quantum phase diagram for t < and V 3> |t| B 
1211 . In terms of the original model, this implies that 5 = 
1/2 dimers with frustrated inter-dimer couplings, J 2 > Ji 
provides an alternative realization of a spin-SS on different 
frustrated lattices. 

As a final example, we consider 5 = 1 Heisenberg model 
with large easy-plane single-ion anisotropy (A = 1 and 
D < in Eq.(l)). For \D\ ^> J the low-energy subspace con- 
sists of the Sf = ±1 states (see Fig.(3jc)). The low-energy 
effective model is once again the t — V Hamiltonian (|4|i with 
t = -J 2 /2D, V = J + J 2 /D and/.* = B- J 2 /D-2n b J, up 
to second order in J/D. rib is the number of bonds per site. 
As in the previous case, a SSphase is realized for V ^> \t\ 
on different frustrated Iatticesyl3,[l2l. The BEC component 
corresponds to spin nematic ordering. 

In conclusion, we have used numerical simulation to study 
the ground state and thermodynamic phase transitions involv- 
ing the spin-supersolid phase in a S = 1 Heisenberg model 
with uniaxial exchange and single-ion anisotropies on a cu- 
bic lattice. The melting of the SS occurs in two steps with 
the XY and Ising ordering disappearing at different temper- 
atures. The transitions are marked by unique features in the 
structure factor and the specific heat which will be useful in 
any experimental detection of the SS. Finally, we discuss sev- 
eral different conditions under which a SS can be realized in 
real spin compunds. 

LANL is supported by US DOE under Contract No. W- 
7405-ENG-36. 
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FIG. 3: (Color online) (a) Square lattice of S=l dimers with an intra- 
dimer Heisenberg AFM interaction Jo and inter-dimer interactions 
Ji and J2. The level diagram shows the low energy subspace of 
the single dimer spectrum in the presence of a magnetic field, (b) 
5 = 1/2 dimers on a triangular lattice and the low energy subspace 
of a single dimer. (c) S = 1 spins on a triangular lattice and the 
energy level splitting for easy-axis single-ion anisotropy. 

For 5=1/2 dimers the low energy subspace of Ho con- 
sists of the singlet and the S z = 1 triplet states in the limit 
Jo <S Ji, J 2 (see Fig. |3lb)). The resulting low-energy effec- 
tive model is a t — V Hamiltonian for hard core bosons: 

Heff - -t Y,( b l b i + b l bi ) + vj2 n wvY. n -> w 

<ij> (iJ> » 

b\ creates a 5 2 = 1 triplet state at site i whereas the singlet 
corresponds to the empty boson state; Uj is the boson number 
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